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1 .  INTRODUCTION 


Let  g  denote  a  complex-valued  function  of  the  real  variable  t,  lntegrable 
with  respect  to  a  Lebesgue  Stleltjes  measure/!  on  an  interval  I.  i5y  the  "quadrature” 
of  g  on  I  we  mean  the  value  of  the  integral 

'l(p)  =  j  g(f)  d/t(f).  (1) 

By "numerical  quadrature"of  g  we  mean  an  (approximate)  evaluation  QM(g)  of  Q(g) 
by  means  of  a  numerical  algorithm  based  on  a  set  of  samples  of  g  on  I  and 
possibly  on  some  of  the  derivatives  of  g  at  the  endpoints  of  I. 

Quadrature  evaluations  occur  in  important  applications.  Some  of  the 
applications  that  we  shall  consider  are: 

(a)  Ordinary  integration  of  g  over  an  interval  I  (/t  is  then  the  ordinary 
Lebesgue  measure  on  the  real  line); 

( b)  Calculation  of  the  nth,  moment  of  g,  where  g  is  a  mass  density,  a 
charge  density,  or  a  probability  density  over  I,  thus 


Q(f)  =  [  (t’)n  c(t')  dt\  (2) 

*  I 

where  n  is  an  arbitrary  nonnegative  integer; 

(c)  Convolution  of  a  signal  g  with  an  impulse  response  h  of  a  dynamical 
system,  thus 


f  h(t-t')  g(t’)  df; 


(3) 


and  (U)  the  Fourier  transform  of  a  signal  g, 


thus 


Q(g)  = 


,to 

e 

oo 


,i2nTt ' 


6(f)  df 


(4) 


The  main  thrust  of  the  work  reported  here  is  the  use  of  "physical  models" 
for  generating  "strictly  numerical"  quadrature  algorithms  as  follows:  We  assume 
that  g  is  generated  by  a  dynamical  system  driven  by  an  imput  u  whose  energy 
is  finite  over  thu  interval  1,  i.o.  ufitT (i).  Through  such  a  selection  of 
a  dynamical  model  it  is  possible  to  incorporate  a-priori  information  about 
g  into  the  design  of  an  algorithm  for  the  numerical  quadrature  of  g. 

Our  model  for  the  generation  of  g  is  of  the  form  (See  Fig. 1 ) : 


F(g, 


u,  '  -  d/dt,  (5) 


where  F  is  a  general  nth.  order  differential  operator.  In  the  linear  case 
we  assume 

*‘(g)  "La  =  (  I  a  Dk)  g=L(D)  g,  J)=d/dt,  (6) 

k=0  R  k=0 

where  a  =  1 ,  and  a,  ,  k  =  0,  ...,  n-1,  are  real  constants.  In  the  nonlinear  case, 

n  k  (n) 

we  assume  that  F  is  a  smooth  function  of  its  (n+1)  arguments  g,  g'  ,  ...,  g  , 

having  continuous  partials  with  respect  to  these  variables  up  to  order  n. 


u  =  F(g)  =  F(S,c\  ....  g(n)) 


Fig.  1  (a)  General  Operator  Model 
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In  both  the  linear  and  nonlinear  cases,  numerical 
quadrature  algorithms  will  be  obtained  as  solutions  of  an  optimization 
problem . 


u 


-jl  Linear  Model 


G 


u  =  p(g)  =  L(D)  r. 


(  Z  a  Dk)  G,  u  =  d/dt 
k-0 


Fig.l(b)  Linear’  Operator  Model 

In  Section  2,  we  consider  the  case  of  a  linear  source  model  for  g  and, 
in  section  we  extend  the  theory  to  the  nonlinear  case. 

Numerical  simulations  performed  at  Rice  show  the  merits  of  the  present 
approach  and  -they  will  be  described  separately. 
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2.  QUADRATURE  BASED  ON  LINEAR  MODELS 

As  we  shall  see  in  subsection  2.1,  linear  models  lead  to  solutions  based  on 
generalized  splines.  Considerable  amount  of  research  on  quadrature  formulas 
based  on  splines  has  been  done  by  Schoenberg  [l]  [2]  ,  Sard  [3]  ,  and  Goloinb  [4] . 

In  what  follows,  we  recast  some  of  this  spline  theory  in  the  language  of 
dynamic  models. 

In  subsections  2.2  through  2.5<  we  show  that  such  models  lead  to  appropriate 
windows  for  use  in  the  numerical  evaluation  of  ordinary  integrals,  moments, 
convolutions,  and  Fourier  transforms.  In  what  we  have  just  stated,  we  have 
assumed  that  the  model  generating  g  is  given  to  us.  If  it  is  unknown ,  we 
show,  in  subsection  2.6,  how  to  retrieve  it  from  the  data  using  a  maximum 
entropy  approach. 

2.1.  Fundamentals  of  Spline  Based  Solutions 
2 

Let  Hn  (l)  denote  the  Sobolev  space  of  complex-valued  functions  g  on  I  =  [a,b] 
such  that  g  ,  k=0 , 1, . . . ,n-l ,  are  absolutely  continuous  and  g^n^£  I) • 

A2/  o 

We  will  denote  by  11^(1)  the  closed,  subspace  of  defined  by 

hn(i)  ~  -CO.  oo):  g(t)  -  0,  t<a,  t  >  b;  g^(a)  =  g^(b)  =  0, 

1.  ...,  n-1;  the  restriction  of  g  to  I 
belongs  to  ir(l)J.  (7) 

Also,  let  L(D)  be  the  linear  operator  defined  by  (6),  i.e. 

n-l  , 

L(D)  =  D  +  £  a.  Dk  .  (8) 

k=0  * 


-  J 


'.la  will  assume  that  the  signal  g  under  consideration  (whose  quadrature 
Is  to  be  obtained)  satisfies  the  following  conditions: 

(i)  g  e  III  (I): 

(ii)  we  are  given  a  set  of  complex  numbers  r^ ,  rM  and  a 


mesh 


a  =  t0<  ll  <  ••• 


<ln  4  Shi  a  b* 


(9) 


with  respect  to  which  g  satisfies  the  interpolating  constraints 

g(  t^)  ~  ^  q  »1—  !»•••»  ^!  (id) 

(iii)  g  belongs  to  the  ellipsoidal  class 


IX">  'L,  ,<r) 

*•  la  (  -00,00  )  J 


(ID 


where  ^ is  some  positive  constant,  and  L(D)  is  assumed  given. 

If  g  satisfies  all  the  above  three  conditions,  we  will  write 

E  t 

Note  that  H  (i)  constitutes  a  Hilbert  space  under  the  inner  product 


-j  (b(D)  g*(t)](l,(D)h(t))  dt,  *  =  complex  conj. 


h2(i) 

n 


We  will  use  the  folllowing  two  wellknown  results  w 


Proposition  1.  The  min -max  problem 


min  max  \\g  -  gjj 


H  (I) 
n  ' 


(1*0 
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has  a  unique  solution  g  which  is  the  L-spllne  defined  by  the  set  of  equations 


and  conditions: 


I.(-D)  J,(D)  e(t)  =  0,  tt  <  t  <  t1+1,  1  =  0 . N;  (15a) 

1  =  1 . X,  (l5b) 

‘•‘C2"'2 (D>  6t«*(D,  (15c) 


where  c"(l)  =  space  of  functions  on  I  with  continuous  derivatives. 


The  u  corresponding  to  such  a  g  is  obtained  simply  by  the  relationship 
u  =  L(I))  g.  /./ 


Proposition  2  .  Let  $  be  a  continuous  linear  functional  on  li  (i). 

n 


The  min-max  problem 


min  max  |$  (g)  -  $>(g)  1 

StifU)  K  t  ^ 


has  a  unique  solution  <£(g)  given  by 

$(C)  =  $(*). 


where  g  is  L-spline  of  Proposition  1. 

Recognizing  that  for  geiljj  (i)  the  quadrature 


(  b 

Q(ft)  =  g(t)  d  ft  (t) 

J  a 


is  a  continuous  linear  functional  on  (l)»  Proposition  2  provides  us  with  a 
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technique  for  calculating  numerical  quadrature  QN(g)  as  a  best  estimate  Q(g) 

of  Q(g)  on  the  basis  of  the  data  g^(a)  =  g^(b)  =  0,  k  =  0,...,  n-1, 

?t( t  )  -  r  i  -  and  the  "physical  model"  that  generates  g,  described 

i  i 

by  the  dynamical  system 

!.([))  g(t)  =  u(t).  (20) 

The  solution  is  to  interpolate  the  data  by  thespline  associated 
with  the  operator  L  and  then  perform  the  quadrature  on  this  spline. 

By  way  of  example,  suppose  we  model  g  as  the  output  of  a  low-pass 
filter  with  half -power  frequency  equal  to  one  (See  Fig.  2).  We  then  have 

L(D)  =  D  +  1,  (21) 

and  the  roots  of  the  characteristic  equations  associated  with  L(D)  and 
L(-D)  are  respectively  -1  and  1. 

Model 


Fig.  2.  Network  Model  for  Basing  Optimal  Numerical 
Quadrature  Algorithm 


-  0  - 


So  we  haves 

g(t>  a+i  e't  +  a\  et  .  ‘i<*  <q+1, 

1  =  0,1 . N  ,  (22) 

where  the  constants  a  ^  and  a  ^  are  associated  with  the  homogeneous 
solutions  pertaining  to  L(D)  and  L(-D),  These  constants  are  determined, 
according  to  (15a)  through  (15c),  by  the  requirements  that  g  satisfy  the  in¬ 
terpolating  constraints  and  the  boundary  conditions,  and  that  g  be  continuous 
on  I. 

In  the  general  case  in  which  the  characteristic  equation 

L(s)  =  0  (23) 

associated  with  the  operator  L(D)  has  n  distinct  roots  s^...,  s  ,  (22) 
generalizes  to 


g(t)  = 


2  [a+ 

j=l 


ij 


sjt 

e  J  + 


ij 


■V] 


ti< 


t  <  t 


i+1, 


i  0,  1,  . . . ,  N ,  (24) 


where  again  the  constants  a  ^  and  a  ^  are  determined  by  the  requirement 

that  the  interpolating  and  boundary  constraints  be  satisfied  in  addition  to  the 
requirement  that  the  derivatives  of  g  up  to  the  order  2n-2  be  continuous  at 
the  knots. 

If  (23)  has  a  repeated  root  say  s  of  multiplicity  ■  V  ,  the  expression 

K  K 

inside  the  square  brackets  in  (24)  for  that  particular  root  is  replaced  by 
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Vk  "  1  s  t 

[2  (aikp  tP  eJ  +  aikp  ^  (25) 

p=0 

,  +  — 

where  a^p  and  are  suitable  constants  that  depend  on  the  data  and 

other  conditions  as  mentioned  above. 

A  set  of  basis  functions  which  makes  the  dependence  of  g  on  the  data 

transparent  is  the  so-called  "cardinal  spline  basis"  [5)  (also  called  "fundamental 

spline  basis"  [6]  )•  This  basis,  which  we  denote  by  B,  is  the  set  of  splines 

t’j(d),...,  b^(t)  satisfying  15(a)  through  15(c)  with  the  data  satisfying 

r«  =  5  . .  (=kronecker  delta)  for  the  spline  b..  Thus  the  spline  b.  has 
x  J  J  J 

zero  values  at  all  the  knots  except  at  the  jth.  knot  where  it  has  value  of 

unity.  It  vanishes  outside  [a,b]  .  The  sketches  of  Fig.  3  exhibit  the 

shapes  of  these  basis  functions.  They  look  and  act  very  much  like  sine  functions 

and  they  apply  to  a  finite  set  of  samples  (rather  than  an  infinite  set  like  the 

sine  functions). 


Fig.  3.  Sketches  of  Cardinal  Spline  Basis  Functions 


-lo¬ 


in  terms  of  the  basis  elements  of  B,  ^  can  be  represented  as 
N 

g(t)  =  ZL  s(ti)  bi(t) 

i=l 

N 

=  Z  ri  bi (t)  •  <26> 

i=l 


We  will  use  this  representation  of  g  in  the  following  subsections. 

Summarizing  the  developments  in  this  subsection,  we  claim  that  we  have 
given  all  the  details  required  to  represent  and  calculate  the  entities  needed 
for  the  quadrature  problems  described  below.  They  were  used  in  the  numerical 
simulations  performed  on  our  computers. 


2.2. 


Ordinary  Integration 

The  first  application  that  we  shall  consider  is  the  numerical  integration 


Q(g) 


g(t)  dt, 


(27) 


't 


1 


based  on  samples  of  g  given  by  (10)  and  on  the  other  conditions  stated  previously. 

Representing  our  optimal  estimate  g  of  g  by  means  of  the  cardinal  spline 
basis  functions  given  by  (26),  we  obtain  the  formula  for  the  "model  based" 
numerical  integration  corresponding  to  (27) ! 

N  N 

QN(s)  =  Z  wi  §(ti)  =  Z  wi  ri  »  (28a) 

i=l  i=l 


where 


-il- 


b^(t)  dt 


i  =1 


,N. 


(28b) 


2.3.  Evaluation  of  Moments 

V/e  remind  the  reader  that  in  this  and  in  the  following  subsections, 
we  assume  the  signal  g  to  vanish  outside  the  interval  (a,b).  This  interval 
is  of  course  allowed  to  be  of  infinite  length  if  necessary,  by  setting 
a  =  -00  and/or  b  =  00. 

An  optimal  algorithm  for  the  numerical  evaluation  of  the  kth.  moment 

i 

now  follows  from  the  preceding  theory  to  be 
N  N 

QN(g)  =  X  wi  s(-t1)  =  H  wi  ri»  (29a) 

i=l  i=l 

where 

tk  b  (t)  dt,  i  =  1 . N.  (29b) 

I 


2.4.  Convolution 

Our  model-based  optimal  numerical  scheme  for  convolution  of  g  with 
an  impulse  response  h  is 


QN(h®g)  (t) 


N 


£ 


"lt  s(V 


N 


X 


w 


it  i 


j 


(30a) 
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h(t-tJ )  b,(f)  dt',  i  =  1 . N. 


(30b) 


2.5.  Fourier  Transforms  and  Filter  Transfer  Functions 

In  the  same  fashion  as  in  the  preceding  subsections,  we  derive  our 
model-based  algorithm  for  the  calculation  of  the  Fourier  transform  of  g  by 


ysH-)  «  z:  “iuJ  g(t1)  = 


N 

i.  w 

i=i 


iwr  i ' 


(31a) 


where 


V,  ■  J”, 


(t)  e'JlJt  dt,  i  =  1 . N. 


(31b) 


Note  that  in  the  special  case  in  which  the  sampling  instants  are  equi¬ 
distant,  if  we  set 

b^t)  =  <f  (t  -  t±), 

(31a)  (  if  we  consider  appropriate  discrete  values  of  oo  )  becomes  the 
discrete  Fourier  transform  of  the  data.  Thus  we  see  that  the  algorithm 
that  we  proposed  is  a  flexible  "mod el -based"  extension  of  the  DFT  algorithm. 

In  the  case  in  which  g  is  the  impulse  response  of  a  filter,  we  conclude 
that  (31a)  and  (31b)  constitute  an  optimal  algorithm  for  the  derivation  of 
its  transfer  function  in  terms  of  the  samples  of  the  impulse  response. 


2.6.  Model  Identification  by  a  Maximum  Entropy  Approach 

Thus  far ,  the  differential  operator  L(D)  describing  the  model  was 
assumed  known.  In  general,  this  may  not  be  the  case,  and  one  would  like 
to  develop  a  methodology  for  determining  L(D)  from  the  data  and  other 
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prior  Information.  Precisely,  this  may  be  formulated  as  follows: 

Problem,  Find  the  constants  a0,a^,..,,a^  ^  in  the  differential  operator 
(8)  from  the  data  r-^  » . . .  ,r^. 

The  identification  procedure  that  we  propose  applies  to  uniformly  sampled 
data  and  is  as  follows.  We  assume  that  the  samples  result  from  an  auto¬ 

regressive  process. 


k+n 


ln-lrk+n-l  +-,,+alrk+l  +  a0rk  ~  uk,  k  =  1 . N_n’  ^ 


where  u^  is  a  white  noise  process.  The  rationale  for  this  assumption  is 
that  it  constitutes  a  discrete  stochastic  version  of  equation  (6). 

From  . .r^,  we  compute  n  samples  of  the  sample  autoavariance  R  (£ ) , 

i  ® 

=  0,.,,n-l.  From  these,  the  coefficients  a.  may  be  obtained  by  requiring  that 

J 

the  Yule-Walker  equations  [7}  be  satisfied  by  R  ,  i.e. 

5 


n-1 

Z  Kg(J-e)  =  -Rg(n-t).  6  =  1 . n  .  (34) 

j=0 

The  coefficients  {a.V  can  also  be  obtained  for  the  process  (33)  using  the 
maximum  entropy  method  in  a  standard  way  [8] . 


3.  QUADRATURE  BASED  ON  NONLINEAR  MODELS 


We  now  propose  a  generalization  of  the  preceding  theory  for  the  case  in 
which  the  signal  g  is  assumed  to  be  generated  by  the  nonlinear  model  (5). 

We  assume  that  g  satisfies  conditions  (i)  and  (il)  stated  at  the  top  of 
page  5  and  (instead  of  belonging  to  the  ellipsoidal  class  (ll) )  belongs  to  the 
class 


(I):  l|F(g,g‘  ...,g(n))||  .  *Y}  (35) 

v  L  (-00, 00)  J  ’ 

where  the  F  is  such  that  Vis  a  convex  and  symmetric  set  in  H2  (I)  and,  in 


addition,  all  nth  order  partials  of  F  with  respect  to  all  its  arguments  are 
continuous. 

jt 

We  recall  from  the  literature  [4]  [9]  that  in  the  linear  case  the  solution 
g  of  (14)  is  also  the  solution  of  the  minimum  norm  problem 


min  (1  L(D)gH2L2(l) 

seH2(i) 

g(ti)=ri'i=1 . N 


(36) 


By  analogy ,  we  seek  the  best  estimate  g  of  g  in  the  nonlinear  case  as  the 
solution  of 


min  !!  F(g)|(  l2(i)  , 

s(ti)=ri»i=1,. . .N 


(37) 
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where  we  recall  that 

F(g)  =  F (g,g* ,...g^)  . 


(38) 


For  simplicity  in  presentation  let  us  assume  that  the  underlying  linear 
spaces  are  over  the  field  of  reals. 

The  solution  g  of  (37)  will  be  called  a  "nonlinear  spline".  Such  g  must 
satisfy  po): 


5 


fZ(g.g' . g(n))  dt  =  0 


(39b) 


^(tf)  =  r^  i=l , . . .  ,N  . 


(39b) 


Performing  the  variation  indicated  in  (39a)  (under  the  fixed  boundary 
conditions  assumed  earlier),  we  are  led  to  the  Euler  equation  (we  omit  the 
details  of  the  derivation) i 
n 


Z(“4t)k  [F(e’Sl 


(40) 


Thus  the  nonlinear  spline  g  is  described  by  the  nonlinear  differential 
equation  (40)  in  each  of  the  subintervals t,  <;  t  <  t  i  =  0,1,...,  N,  and  by. 

i+i , 

the  conditions  (15b)  and  (15c). 

The  optimal  numerical  quadrature  based  on  the  nonlinear  model  developed 
above  is  simply  given  by  (19)  with  g  replaced  by  g  of  the  preceeding  paragraph. 


4.  CONCLUSION 


A  physical  modeling  basis  and  framework  have  been  provided  for  constructing 
various  quadrature algorithms .  For  this  purpose,  both  existing  and  new  results 
on  spline  theory  have  been  presented. 
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